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Part  of  the  Project  RAND  research  program  consists  of 
basic  supporting  studies  in  mathematics.  In  this  Memorandum 
the  authors  discuss  a  method  for  solving  large  systems  of 
differential  equations  where  the  solution  is  subject  to 
certain  boundary  conditions. 


SUMMARY 


Two— point  boundary— value  problems  for  second— order 
systems  of  linear  differential  equations  are  usually 
solved  by  a  process  involving  the  inversion  of  a  certain 
matrix.  If  the  system  is  too  large ,  it  may  be  difficult 
to  compute  this  inverse  to  a  high  degree  of  accuracy. 

The  purpose  of  this  paper  is  to  demonstrate  that  this 
difficulty  can  in  some  cases  be  circumvented  by  applying 
a  method  like  that  of  Bodewig  and  Hotelling. 
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ON  THE  COMPUTATIONAL  SOLUTION  OP 
TWO-POINT  BOUNDARY-VALUE  PROBLEMS 

1 .  INTRODUCTION 

Consider  (as  in  [ li  )  the  n— dimensional  vector 
differential  equation 

(1.1)  x"  +  A(t)x  =  0  , 

where  the  solution  is  subject  to  the  boundary  conditions 

(1.2)  x(0)  =  c,  x(l)  =  d. 

The  problem  is  generally  solved  as  follows.  Let 
and  X2  denote  the  matrix  solution  of 

(1.3)  X"  +  A(t)X  =  0 
satisfying  the  initial  conditions 

(1.4)  X1(0)  =  I,  X'(0)  =  0, 

Xg (0)  =  0,  X£(0)  =  I. 

If  g  represents  the  (unknown)  value  of  x'(0), 
where  x(t)  is  the  solution  to  the  problem,  then 

(1.5)  g  =  X2(l)“1[d  -  X1(l)c) . 

If  Xg(l)  is  singular,  then  there  may  be  many 
solutions,  or  none,  and  (1.5)*  of  course,  makes  no  sense. 

If  n  is  large,  it  may  be  difficult  to  compute 
X^(l)  to  a  high  degree  of  accuracy.  The  purpose  of  this 
paper  is  to  discuss  a  method  of  overcoming  this  difficulty. 


2.  AN  ITERATIVE  TECHNIQUE 


Let  X2(l)  be  some  approximation  to  X2*(l), 

Define 

(2.1)  &1  =  XgdJtd  ~  X1(l)c], 

gn  =  X*(l)[d  -  X1  (l)c  -  X2(l)grv_1]  +  grv_r 
Then  we  have  the  following  theorem: 

Theorem.  If  the  spectral  radius  of  I  -  X2(l)Xp(l) 

is  less  than  one,  then  the  sequence  (g„)  defined  by 
- — - * — • —  n  - 

(2.1)  converges  to  g,  the  unique  solution  of  (1.5) • 

Proof .  First  note  that  if  I  -  Xg  ( 1 )  Xg  ( 1 )  has 
spectral  radius  less  than  one,  then  Xg(l)X2(l)  must  be 
nonsingular.  Thus  X2(l)  and  X2(l)  are  nonsingular, 
which  means  that  (1.5)  has  a  unique  solution.  If  g  is 
the  unique  solution  of  (l.5)>  then 

(2.2)  gn  -  g  =  Xg ( 1 ) [  d  -  X1(l)c  -  X2(l)gr>_1]  +  gn_x  - 

=  X2(l)[d  -  X1(l)c  -  X2(l)gn_]_] 

-  Xg(l)Id  -  X1(l)c  -  X2(l)g]  +  -  g 

=  (i  -  xj(i)x2(i))(gn_1  -  g). 

*X* 

If  the  spectral  radius  of  I  -  X2(l)X2(l)  is  less 
than  one,  this  shows  that  (gn  -  g)  goes  to  zero  as  n 
goes  to  infinity,  and  this  concludes  the  proof.  This 


-3- 


theorem  may  be  viewed  as  an  application  of  a  method  of 
matrix  inversion  like  that  of  Bodewig  and  Hotelling  (see 
[3],  [4]  for  additional  references). 

2 

Corollary.  If  A(t)  =  B  j  a  constant  positive- 
definite  matrix,  then  taking  X^(l)  =  Xg(l)  makes  (gn) 
converge  to  the  solution. 

Proof.  Since  X2(l)  =  B-'1'  sin  3,  it  follows  that 
the  eigenvalues  of  Xg(l)  all  have  absolute  value  less 

O 

than  onej  and  thus  all  the  eigenvalues  of  X2(l)  are 
between  0  and  one. 

Corollary.  If  each  element  of  I  -  X2(l)X2(l)  is 
less  in  absolute  value  than  l/n,  then  (g  }  converges 
to  the  solution. 

2 

Corollary.  If  A(t)  =  -  B  ,  where  ‘B  is  a  matrix 
with  only  real  eigenvalues  each  of  which  is  greater  than 
zero,  then  taking  X2(l)  =  2Be~B  makes  (g^j  converge 
to  the  solution. 

Bt  -Bt 

Proof.  X2(t)  =  B  ( - 2—^ - ),  whence 

Xp(l)X2(l)  equals  I  -  e“2j3. 

Corollary.  If  ( t ) ,  Y0 ( t )  are  solutions  to 
Y"  +  A(1  —  t)Y  =  0  satisfying  initial  conditions  like 
(1.4),  then  taking  X2(l)  =  Yj(i)  will  make  (gn) 
converge  to  the  solution  if  Y£  ( 1 ) XI,  ( 1 )  has  spectral 


radius  less  than  one. 
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Proof.  Y£(l)X£(l)  =  I  -  Yj(l)X2(l). 

Corollary.  If  Xgtl)  =  dA,  where  A  Is  the 
transpose  of  X2(l)  and  d  Is  a  positive  constant 
chosen  to  be  less  than  twice  the  reciprocal  of  the 
sum  of  the  absolute  values  of  each  row  of  A Xg ( 1 ) ,  then 
( g.  }  converges  to  the  solution . 

Mote  that  this  last  corollary  is  not  apt  to  be 
computationally  useful,  however,  since  if  X2(l)  has 
some  very  small  eigenvalues  (and  thus  is  hard  to  invert), 
under  the  above  procedure  I  -  X2(l)X0(l)  will  have 
spectral  radius  very  close  to  one,  so  that  convergence 
will  be  slow. 
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